Microscopic nonlinear optical response: Analysis and calculations with the Floquet–Bloch formalism

We analyze microscopic nonlinear optical response of periodic structures within the Floquet–Bloch formalism. The analysis is focused on the real-space distributions of optically induced charge and electron current density within the unit cell of a crystal. We demonstrate that the time-reversal symmetry of a crystal determines the phases of the temporal oscillations of these distributions. We further analyze their spatial symmetries and connection to macroscopic optical response. We illustrate our study with ab initio calculations that combine density functional theory with the Floquet–Bloch formalism. The calculations provide time-dependent optically induced charge distributions and electron current densities within the unit cells of a crystal with inversion symmetry MgO and a crystal without inversion symmetry GaAs in response to a strong-field excitation. The real-space, microscopic view on nonlinear optical response provides insightful information about the strong field–matter interaction.


I. INTRODUCTION
Strong-field excitation by light can be used to induce various important mechanisms in solids, such as manipulation of electronic gaps and structure by light, [1][2][3][4][5][6][7][8][9][10][11][12] or generation of high harmonics (HHG). 13Such processes, on the one hand, have a big potential for the development of petahertz electronics and, on the other hand, raise many scientific questions about the mechanisms behind them. 14,15][18][19][20][21][22] The optical response of crystals has been extensively investigated for more than a hundred years.Such studies have predominantly concentrated on the macroscopic optical response of a crystal, since it determines typical experimentally detectable observables, such as harmonic generation.The macroscopic optical response of a crystal in most cases results from an induced dipole moment.Typically, the radiation power produced by the oscillating dipole moment dominates over the radiation power produced by higher-order moments in the optical regime, 23 and the macroscopic polarization in a dielectric material results from the induced dipole moment. 24For this reason, it is customary to relate linear and nonlinear optical response to the induction of dipole moments.
In recent years, laser-driven electron dynamics has gained considerable attention due to remarkable achievements in the field of nonlinear optics including the generation of high harmonics (HHG) in solids, 13,25 optical-field-induced currents in dielectrics, 1 manipulation of electric properties of a dielectric with the electric field of light, 2 control of coherent Bloch oscillations, 5,6 subcycle terahertz nonlinear optical effects, 3 and coherent control of currents in semiconductors using synthesized optical waveforms. 426][27][28][29][30][31][32][33][34][35][36][37][38] We employ the Floquet-Bloch framework, which provides us with a material-specific description of light-matter interaction beyond perturbation theory, [39][40][41] and is also a convenient tool to analyze properties of laser-dressed systems, such as selection rules of HHG in periodic structures. 30,39,42,43Here, employing the Floquet-Bloch formalism, we analyze spatial symmetry and temporal behavior of microscopic optically induced charge and electron-current distributions.
The Floquet formalism implies that the electric field of the driving field is temporally periodic.It has been shown in Refs.44 and 45 that this approximation is already justified for a strong-field optical field with a duration comprising several tens of optical cycles.Throughout this paper, we consider a laser-dressed crystal in a state that is characterized by a single Floquet state.To justify this approximation, it should be assured that the driving field is not too strong to bring the laser-dressed system into a superposition of several Floquet states. 46he applicability of these approximations can be verified with the radiation spectrum generated by the system through the driving field.If they do not hold, the radiation spectrum will have additional peaks besides the harmonic ones.Thus, our study applies to the regime in which each generated radiation peak can be clearly assigned to an integer multiple of the driving-laser frequency.
The article is organized as follows.We analyze in Sec.II temporal and spatial properties of the microscopic nonlinear optical response within the Floquet-Bloch formalism.Namely, we focus on analysis of optically induced charge distribution and electron current density in real space.In Sec.III, we provide details of our calculation scheme based on the density functional theory (DFT) combined with the Floquet-Bloch formalism to calculate microscopic nonlinear optical response.We calculate the microscopic optical response of two bandgap materials, the insulator MgO, and the semiconductor GaAs driven by an intense optical field.By studying these two prototypical materials, we cover two types of crystals, one without inversion symmetry, GaAs, and one with inversion symmetry, MgO.

II. MICROSCOPIC OPTICALLY INDUCED CHARGE AND ELECTRON-CURRENT DISTRIBUTIONS
47][48][49] The Hamiltonian of a laser-dressed crystal is given by (1) Here, Ĥ el is the mean-field Hamiltonian of the unperturbed crystal with one-body eigenstates ju mkr i, where k is the Bloch wave vector, m is the band, and r is the spin index.V c ðrÞ ¼ V c ðr þ RÞ is a spaceperiodic crystal field potential, and R is a lattice vector.According to the Bloch theorem, 50 the corresponding one-body wave function of ju mkr i has the form u mkr ðrÞ ¼ e ikÁr u mkr ðrÞ, where u mkr ðrÞ ¼ u mkr ðr þ RÞ is a space-periodic function.p is the canonical momentum of an electron, and ŵ † ( ŵ) is the electron creation (annihilation) field operator.Ĥ em is the Hamiltonian of the electromagnetic field, and Ĥ int the interaction Hamiltonian between the electromagnetic field and the electronic system, which we describe within the dipole approximation.â † j;s (â j;s ) creates (annihilates) a photon with wave vector j and polarization s.We assume that only the j 0 , s 0 mode with a corresponding polarization vector 0 and the energy x ¼ jj 0 jc, where c is the speed of light, is occupied in the driving electromagnetic field, and that the state of the field is described by a single-mode coherent state ja; ti.Âem ðrÞ is the vector potential operator of the electromagnetic field, and a is the fine-structure constant.We use atomic units for these and the following expressions.
Since the state of the electromagnetic field ja; ti is unaffected by the interaction with the electronic system by assumption, the one-body solution of the time-dependent Schr€ odinger equation idjw i;kr ; ti=dt ¼ Ĥ elÀem jw i;kr ; ti can be represented as jw i;kr ; ti ¼ j/ el i;kr ; tija; ti with the corresponding electronic one-body Floquet-Bloch wave function 39,40,48 / el i;kr ðr; tÞ ¼ X m;l c i m;kr;l e Àilxt u mkr ðrÞ; where l is an integer and the c i m;kr;l are expansion coefficients.

A. Optically induced charge distributions
Within the Floquet formalism, the electron density of the laserdressed system evolves in time as 49 qðr; tÞ ¼ X l e ilxt ql ðrÞ; with lth-order optically induced charge distributions c iÃ m 0 ;kr;l 0 þl c i m;kr;l 0 u † m 0 kr ðrÞu mkr ðrÞ; (7)   where V uc is the volume of a unit cell, i denotes the index of occupied one-body Floquet-Bloch states, and the integration is over the Brillouin zone.The lth-order optically induced charge distributions can be connected to properties of the lth-order macroscopic optical response.For example, the polarization is determined by the dipole moment of the time-dependent electron density Thus, the amplitudes ql ðrÞ are optically induced charge distributions that give rise to a lth-order macroscopic optical response.
Our connection of the optically induced charge distributions derived within the Floquet-Bloch formalism to the macroscopic polarization is consistent with the perturbative derivation of macroscopic polarization by Bloembergen and Shen. 51,52Their expansion of the polarization in orders of x, P ¼ Pð1Þ ðxÞ þ Pð2Þ ð2xÞ þ Pð3Þ ð3xÞ þ Á Á Á also holds within the Floquet formalism.However, lth-order components derived within the Floquet formalism do not have to scale as the lth power of the electric-field amplitude, in contrast to the ones derived within the perturbation theory by Bloembergen and Shen.
The volume integral of a lth-order optically induced charge distribution over a unit cell is given by ð d 3 rq l ðrÞ ¼ where the sum goes over spin directions r.We now apply the orthogonality of the expansion coefficients P m;l 0 c iÃ m;kr;l 0 þl c i m;kr;l 0 ¼ P m;l 0 d l 0 þl;l 0 and obtain ð d 3 rq l ðrÞ ¼ N el d l;0 : The volume integral of the optically induced charge distributions of a nonzero order that enter in the time-dependent part of qðr; tÞ is zero, Ð d 3 rq l6 ¼0 ðrÞ ¼ 0. These relations indicate that the interaction of a crystal with light leads to a dynamical redistribution of charges, which has positive and negative regions relative to the field-free electron density.These positive and negative regions coherently oscillate and cancel each other when volume integrated.The positively charged regions are due to electron holes in valence bands, and the negatively charged regions are due to electrons in conduction bands.

B. Time-reversal symmetry and Floquet-Bloch wave functions
In the following, we will show that the time-reversal symmetry of a crystal determines the temporal behavior of the laser-driven spinaveraged electron oscillations.The derivations below are applied to a general situation in which the crystal is not necessarily invariant under inversion symmetry in real space.First, let us consider how time-reversal symmetry influences the properties of the Floquet-Bloch functions.For a Bloch function u mk" ðrÞ of an unperturbed crystal that obeys timereversal symmetry, it is valid that u mk" ðrÞ ¼ u Ã mÀk# ðrÞ. 50Similar to the proof of this property, we show below that a one-body electronic wave function of a laser-dressed crystal with time-reversal symmetry obeys / el i;k" ðr; tÞ ¼ / Ãel i;Àk# ðr; T=2 À tÞ; where T ¼ 2p=x is the period of the driving electromagnetic field.
Applying to the time-dependent Schr€ odinger equation that the coherent state ja; ti obeys idja; ti=dt ¼ Ĥ em ja; ti, and multiplying Eq. ( 13) by ha; tj, we obtain The matrix element ha; tj Ĥ int ja; ti gives the interaction Hamiltonian in the classical limit where A em ðr 0 ; tÞ ¼ ðc=xÞE em ðr 0 Þ cosðxtÞ with E em ðr 0 Þ being the amplitude of the electric field.Since we apply the dipole approximation, we ignored the spatial variation of the vector potential and the electric-field amplitude, and substituted the position of the crystal r 0 for r in A em ðr 0 ; tÞ and E em ðr 0 Þ.Thus, the one-body electronic wave function obeys In order to prove Eq. ( 12), we take the complex conjugate of Eq. ( 15) and use that We apply the operator rs ¼ ir y , where ry is a Pauli matrix, which acts only on spin on both sides of the equation above We further use that Ĥ cl int ðtÞ commutes with rs , and Ĥ Ã el ¼ r † s Ĥ el rs for crystals with time-reversal symmetry, 50 which results in We now rewrite the above expression for the time T=2 À t taking into account that A em ðr 0 ; T=2 À tÞ ¼ ÀA em ðr 0 ; tÞ and that the derivative over time changes it sign for a negative time Thus, we obtain that if the wave function / el i;k" ðr; tÞ is a solution of the time-dependent Schr€ odinger equation, then / elÃ i;k# ðr; T=2 À tÞ is also a solution.Since / el i;k" ðr þ R; tÞ ¼ e ikÁR / el i;k" ðr; tÞ, it must be valid that rs / elÃ i;k" ðr þ R; T=2 À tÞ ¼ e ÀikÁR rs / elÃ i;k" ðr; T=2 À tÞ.Thus, the solution / elÃ i;k# ðr; T=2 À tÞ is the solution at the Bloch vector Àk, and we have proven Eq. ( 12).Equation ( 12) leads to a connection between the corresponding expansion coefficients of the Floquet-Bloch functions [cf.Eq. ( 5)] which follows from the phase relation e ilxt ¼ ðÀ1Þ l e ÀilxðT=2ÀtÞ .The relation in Eq. ( 20) is necessary for the analysis of the temporal dependence of optically induced charge distributions and electron current densities.

C. Oscillations of optically induced charge distributions
We now make use of Eq. ( 20) and the property of the Bloch functions of crystals with time-reversal symmetry that u mk" ðrÞ ¼ u Ã mÀk# ðrÞ 50 to connect terms with opposite k in the integral for ql ðrÞ in Eq. ( 7).This allows us to reduce the integration over the Brillouin zone to half of the Brillouin zone (HBZ) for the calculation of spin-averaged electron density in the absence of spin-orbit coupling (SOC), which lifts the degeneracy of the spin components of electronic states It follows from this relation that even-order optically induced spinaveraged charge distributions ql even ðrÞ are real functions, whereas oddorder optically induced spin-averaged charge distributions ql odd ðrÞ are purely imaginary.Spin-resolved charge distributions in the presence of spin-orbit coupling have both nonzero real and imaginary parts, which are connected via ql odd " ðrÞ ¼ Àq Ã l odd # ðrÞ and ql even " ðrÞ ¼ qÃ l even # ðrÞ.However, upon spin averaging, the imaginary or real parts of even or odd amplitudes, respectively, vanish.In the following, electron densities and electron current densities are spin-averaged unless stated otherwise.
The time-reversal symmetry has an important consequence for the time dependence of the electron density.We use Eq. ( 21) to combine terms with opposite l in the expression for the time-dependent electron density in Eq. ( 6).It can also be easily derived that ql ðrÞ ¼ qÃ Àl ðrÞ independently of the crystal symmetry.This results in the time dependence of the electron density qðr; tÞ ¼ q0 ðrÞ À . 1 ðrÞ sinðxtÞ þ . 2 ðrÞ cosð2xtÞ À .
. l even cosðl even xtÞ; for the driving field with the vector potential evolving as cosðxtÞ.In Eq. ( 22), we redefined the optically induced charge distributions as follows: . l even ðrÞ ¼ 2Re ql even ðrÞ Â Ã ; . l odd ðrÞ ¼ 2Im ql odd ðrÞ Â Ã ; which are real functions for both even and odd l.This representation of optically induced charge distributions is more insightful in comparison with the functions ql ðrÞ, since it demonstrates the actual time dependence of the light-induced charge distributions.We, thus, obtain that time-reversal symmetry determines the phases of lth-order oscillations of optically induced charge distributions.Broken time-reversal symmetry would lead to a different relation between the phases of induced-charge and electric-field oscillations.

D. Oscillations of electron-current-density amplitudes
We now analyze time dependence of the electron current density.The electron current density in the presence of the electromagnetic field is given by 39,53 jðr; tÞ ¼ Àqðr; tÞA em ðr 0 ; tÞ / elÃ i;kr ðr; tÞr/ el i;kr ðr; tÞ Using the relation between the expansion coefficients of the Floquet-Bloch functions at opposite k due to the time-reversal symmetry in Eq. (20), we obtain the time evolution of the electron current density jðr; tÞ ¼ À X l odd !1 j l odd ðrÞ cosðl odd xtÞ À X l even !2j l even ðrÞ sinðl even xtÞ; where with are real-valued amplitudes of the electron current density.Thus, the electron-current-density amplitudes oscillate with a phase shifted by p=2 with respect to the oscillations of the charge distributions of the same order.When the absolute value of the lth-order charge oscillation is at a maximum, the lth-order electron current density is zero and vice versa.
The expression of the electron current density in Eq. ( 26) is presented in velocity gauge.As shown in Ref. 54, the electron current density calculated in velocity gauge equals the electron current density in length gauge taking the sum of the interband and intraband contributions into account.The separate contributions of interband and intraband currents, however, must be corrected to be consistent in both gauges.Here, we discuss the total electron current density, not separate contributions.
Applying the continuity equation divjðr; tÞ ¼ À@qðr; tÞ=@t l even x.l even ðrÞ sinðl even xtÞ; we obtain the following relation between the amplitudes of the electron current density and the amplitudes of the electron density: Let us now analyze the connection between the dipole moment of the optically induced charge distributions Ð d 3 rr.l ðrÞ and the electron-current-density amplitudes ð Equation ( 33) follows from vector-algebra relations for the dipole moment of a divergence.If the electron current density j l ðrÞ on the boundary of a unit cell is zero, the first term in Eq. ( 33) disappears, and the lth-order macroscopic polarization is proportional to the volume integral of the lth-order electron-current-density amplitudes

E. Optically induced charge distributions and electron current density in a crystal with inversion symmetry
In this subsection, we assume that the crystal exposed to periodic driving is invariant under inversion symmetry and, consequently, Ĥ el ðrÞ ¼ Ĥ el ðÀrÞ, and derive the symmetry properties of electron density and electron-current-density amplitudes.The interaction of Hamiltonian Ĥ cl int is invariant under the transformations r !Àr and t !t À T=2.Thus, it follows from the time-dependent Schr€ odinger equation in Eq. ( 15) that: The Bloch functions of crystals that are invariant under inversion symmetry obey the relation u mkr ðrÞ ¼ u mÀkr ðÀrÞ ¼ u Ã mkÀr ðÀrÞ. 50Substitution of the relations between Bloch and Floquet-Bloch functions into the expansion of the Floquet-Bloch functions in Eq. ( 5) gives a relation between the expansion coefficients at opposite k, Comparing it with the relation between the coefficients due to timereversal symmetry in Eq. ( 20), we find that the coefficients c i m;kr;l are real.
Substitution of these properties into the expression for the optically induced charge distributions via Bloch functions in Eq. ( 7) leads to the following connection between complex amplitudes at opposite r: The property of these amplitudes that either their imaginary or real part is zero depending on the parity of l determines how they behave under inversion symmetry.The same holds for the real-valued representation of the optically induced charge distributions in Eqs. ( 23) and ( 24) that all even-order optically induced charge distributions are invariant under inversion symmetry, whereas all odd-order optically induced charge distributions are opposite under inversion symmetry .l even ðrÞ ¼ .l even ðÀrÞ; . l odd ðrÞ ¼ À. l odd ðÀrÞ: Analogously, using that u mkr ðÀrÞ ¼ u Ã mkÀr ðrÞ, 50 r Àr ¼ Àr r and the coefficients c i m;kr;l being real, we obtain the following symmetry properties of the current-density amplitudes: j l even ðrÞ ¼ Àj l even ðÀrÞ; (41) j l odd ðrÞ ¼ j l odd ðÀrÞ: The volume integral of the even-order current density amplitudes is zero ð d 3 rj l even ðrÞ ¼ 0; (43)   in agreement with the selection rule that even-order harmonics from crystals invariant under inversion symmetry are forbidden. 55

III. MICROSCOPIC OPTICAL RESPONSE IN BANDGAP CRYSTALS MgO AND GaAs A. Computational details
We diagonalize the Floquet-Bloch Hamiltonian as described in Refs.39 and 49.We calculate the one-body wave functions u mkr of the field-free Hamiltonian Ĥ el with the density functional theory using the ABINIT software package [56][57][58] in combination with Troullier-Martins pseudopotentials. 59The functions u mkr of valence bands and conduction bands are calculated on a dense grid of k points in half of the Brillouin zone.The numbers of blocks of the Floquet-Bloch matrix, k points, and bands are increased in the computations until convergence of the Fourier components of the optically induced charge distributions is reached.We apply the scissors approximation 60 to correct a bandgap to its experimental value.
The conduction bands that are necessary to converge the optical response are actually the bands into which electrons are excited by the electromagnetic field with a nonvanishing probability.The number of conduction bands involved in the interaction with the optical field strongly depends on the intensity of the optical field and crystal properties.This number increases with the intensity of the optical field and is well above ten in the nonperturbative regime.There are several reasons for such a high required number of conduction bands.The first reason is that the higher the intensity of the optical field, the larger is the probability of an off-resonant transition into energetically high conduction bands.For example, a transition from a valence band into a conduction band with an energy difference detuned by 10 eV from the photon energy of the driving field can contribute to the first harmonic, if the intensity of the optical field is sufficiently high.
The next reason is that the higher the intensity of the optical field, the larger is the probability of a resonant multiphoton transition.As an example, let us consider the seven-photon absorption process induced by a field with a photon energy of x ¼ 1:55 eV.A transition from a valence band to a conduction band with an energy difference of 7x ¼ 10:85 eV is resonant and should have a dominating contribution to this process.Then, it is crucial to take into account the conduction bands lying at % 11 eV above the outermost valence band for the calculation of the seventh-order optical response.The number of generated harmonics increases with increasing field intensity, and so should increase the number of conduction bands that are necessary to take into account for the calculation of a high harmonic spectrum.Multiphoton transitions also contribute to the optical response at lower orders.For example, seven-photon absorption combined with six-photon emission can contribute to the first-order response.Such contributions would become noticeable under conditions in which perturbation theory is likely to be no longer valid.We used the same computational scheme for the calculation of an x-ray-optical wave-mixing signal in Ref. 49.We have shown that x-ray-optical wave-mixing signal is connected to the Fourier transform of optically induced charge distributions from the real space to the momentum space.That is why in that study, we focused on the calculation of optically induced charge distributions in the momentum space and did not analyze their real-space properties.

B. Crystal with inversion symmetry, MgO
We first calculate the microscopic optical response of a crystal with inversion symmetry, MgO.We consider driving optical field with an intensity of I em ¼ 2 Â 10 12 W/cm 2 , a photon energy of 1.55 eV, and polarization axis ¼ ð0; 0; 1Þ.The calculation is performed using a 24 Â 24 Â 24 Monkhorst-Pack grid, four valence and sixteen conduction bands, and 81 blocks of the Floquet Hamiltonian, which are necessary to reach convergence.The driving field and the computational parameters are the same as we used to calculate the subcycle-unresolved x-ray-optical wave-mixing signal from laserdressed MgO in Ref. 49. There, we showed that an optical field of 2 Â 10 12 W/cm 2 drives electron dynamics in MgO nonperturbatively.For the current computation, we additionally apply the scissors approximation 60 to correct the bandgap from the calculated 5.6 eV to the experimental value of 7.8 eV. 61igure 1 shows the calculated microscopic response of the MgO crystal depending on the phase of the driving electromagnetic field with the electric field evolving as E em sinðxtÞ.A cut of a unit cell centered at the Mg atom is shown.The first column displays the firstorder oscillations of the electronic state, i.e., the oscillations of the electron density and the electron current density with frequency x in response to the driving electromagnetic field.As shown in Sec.II B, the first-order oscillations of the electronic state of the laser-driven crystal comprise the oscillations of the electron current density as Àj 1 ðrÞ cosðxtÞ and the oscillations of the electron density as À. 1 sinðxtÞ.
The first-order optically induced electronic state at xt ¼ 0 is shown in Fig. 1(a).At this phase, the electric field of the optical field is zero, the magnitude of the microscopic first-order electron current is at a maximum, and the first-order electron density is zero.Thus, Fig. 1(a) shows only the first-order electron current density Àj 1 ðrÞ.The lth-order electron current amplitudes j l ðrÞ are three-dimensional vector fields that are nonzero at most points within the unit cell.For the purpose of intuitive visual representation, the electron current densities in this and other figures are plotted on a sparse grid and only vectors with magnitudes jj l ðrÞj larger than a certain minimum threshold are shown.The magnitudes of j l ðrÞ are color-coded, and their values are in atomic units.The minimum threshold for jj l ðrÞj in a given plot is the minimum value of the corresponding color box.
The first-order electron current density in Fig. 1(a) points along the driving-field polarization direction, and its magnitude jj l ðrÞj has pronounced peaks at the oxygen atoms.The electron current causes the charge to rearrange within the unit cell and, at xt ¼ p=2, the firstorder electron density À. 1 sinðxtÞ reaches the maximal magnitude.
Figure 1(b) shows the first-order electron density at xt ¼ p=2, when the electric field is at a maximum and the first-order electron current density is zero.The electron densities in this and other figures are represented in terms of an isosurface using VESTA software. 62The yellow and blue colors, respectively, represent negative and positive charges relative to the field-free electron density.As shown in Sec.II E, the optically induced charge distributions of a laser-driven crystal with inversion symmetry also have inversion symmetry.Thereby, oddorder optically induced charge distributions are antisymmetric with respect to a center of symmetry.Consistently, the values of the firstorder electron density in Fig. 1(b) at positions r and Àr relative to the Mg atom are opposite.The positive and negative charges in Fig. 1(b) alternate along the z axis parallel to the optical-field polarization.This charge alignment is consistent with the macroscopic first-order polarization of MgO being aligned with the electric field as expected.
At xt ¼ p, the magnitude of the first-order electron current density is again at a maximum, but opposite to that at xt ¼ 0. The electron current now causes the charge to redistribute in the opposite direction.At xt ¼ 3p=2, the first-order electron density is opposite to the first-order electron density at xt ¼ p=2.The macroscopic polarization is again aligned with the electric field.
The second column of Fig. 1 shows the second-order oscillations, i.e., the oscillations with frequency x 2 ¼ 2x, of the electronic state in response to the driving field.The second-order oscillations of the electronic state of the laser-driven crystal are shown, comprising the oscillations of the electron current density as Àj 2 ðrÞ sinð2xtÞ and the oscillations of the electron density as . 2 cosð2xtÞ.The second-order macroscopic polarizability tensor of the MgO crystal is zero because of inversion symmetry.Surprisingly, we find that the second-order microscopic optical response is nonzero.In order to understand this phenomenon, let us look into the second-order electronic state at x 2 t ¼ 0 shown in Fig. 1(e).At this phase, the second-order electron current density is zero and the magnitude of the second-order electron density is maximal.
As discussed in Sec.II E, an even-order optically induced charge distribution of a centrosymmetric crystal is also centrosymmetric.Consistently, the second-order electron density in Fig. 1(e) is centrosymmetric with respect to the Mg atom.Therefore, the charge distribution in Fig. 1(e) does not have a dipole moment and results in zero macroscopic polarization.Analogously, the general statement that an even-order optically induced charge distribution of a centrosymmetric crystal is also centrosymmetric is consistent with the selection rule that the even-order macroscopic polarization of such a crystal is zero.
Figure 1(f) shows the second-order electronic state at xt ¼ p=4, when the magnitude of the second-order electronic current density is at a maximum.Consistently with the discussion in Sec.II E, the microscopic second-order electron current density is antisymmetric with respect to the center of symmetry.
The third column of Fig. 1 shows the third-order oscillations of the electronic state of laser-driven MgO crystal, comprising the oscillations of the electron density as À. 3 sinð3xtÞ and of the electron current density as Àj 3 cosð3xtÞ.Figure 1(i) shows the thirdorder electronic state at xt ¼ 0, which is given by the third-order electron current density.Like the first-order electron current density in Fig. 1(a), j 3 ðrÞ points predominantly in the direction of the driving-field polarization.Figure 1(j) shows the third-order electron density at xt ¼ p=6.Interestingly, it has a very similar structure to the first-order electron density in Fig. 1(b).The charge distribution in Fig. 1(j) also clearly indicates that the third-order macroscopic polarization points along the driving-field polarization direction.
The fourth column of Fig. 1 shows the fourth-order oscillations of the electronic state that comprise the oscillations of the electron density as .4 cosð4xtÞ and of the electron current density as Àj 4 sinð4xtÞ.The fourth-order electron density at xt ¼ 0 shown in Fig. 1(m) has a very similar structure to the second-order electron density in Fig. 1(e).Since the fourth-order charge distribution is centrosymmetric, it leads to zero fourth-order macroscopic polarization.
It is possible to access optically induced charge distributions in an experiment.They can be probed by x-ray-optical wave mixing, in which x-ray diffraction is performed during the action of optical excitation of a crystal. 49,63,64We demonstrated that the x-ray-optical wave-mixing signal encodes the square of Gth Fourier component of lth-order optically induced charge distributions, j Ð d 3 re iGÁr .l ðrÞj 2 . 49They determine the intensity of the lth-order side peaks to a Bragg peak at the reflection G that is the signal separated by lx in energy and by lq in the momentum space from this Bragg peak, where q is the momentum of optical photon.In that study, we presented the calculation of an x-ray-optical wave-mixing signal from a laser-driven MgO as an example.We found several surprising phenomena that we could not explain.We can now understand them looking at the figures showing the microscopic optical response of MgO.
The first unexpected result was that the intensities of even-order side peaks that are proportional to j Ð d 3 re iGÁr .l even ðrÞj 2 were nonzero, although even-order harmonics of MgO are zero.We now found that the microscopic even-order optical response of MgO is nonzero although it results in zero macroscopic optical response.In the x-rayoptical wave mixing experiment, x rays give access to the atomic scale and reveal the microscopic optical response.
The second surprising observation was that the intensities of the odd-order side peaks that are proportional to j Ð d 3 re iGÁr .l odd ðrÞj 2 were zero at G perpendicular to the driving-field polarization , whereas the intensities of the even-order side peaks did not change considerably at G ? .Comparing the odd-order and even-order amplitudes of the electron density, and the odd-order and even-order amplitudes of the electron current density in Fig. 1, this behavior becomes clear.The odd-order electron current densities are aligned along the driving-field polarization direction, such that Ð d 3 re iGÁr .l ðrÞ ¼ Ð d 3 re iGÁr divj l ðrÞ ¼ G Á Ð d 3 re iGÁr j l ðrÞ is zero at G ? .The even-order density amplitudes are close to a spherically symmetric distribution, and their Fourier components do not have such a pronounced angular dependence.That is why their Fourier components do not strongly depend on the angle between G and .

C. Crystal without inversion symmetry: GaAs
We consider the microscopic optical response of a crystal without inversion symmetry, GaAs, which has a bandgap of 1.42 eV. 65We consider optical excitation by an optical field of 1 eV photon energy and an intensity of 4 Â 10 11 W/cm 2 .The factor ffiffiffiffiffiffi I em p =x entering the off diagonal matrix elements of the Floquet Hamiltonian 39 is the same as in the calculation of optical response of MgO.28 Â 28 Â 28 Monkhorst-Pack grid, 4 valence and 56 conduction bands, and 151 blocks of the Floquet-Bloch Hamiltonian are necessary to converge the results.According to the second-order susceptibility tensor of the space group F 43m, 55 the second-order macroscopic response of GaAs driven by a field polarized along the (1, 1, 1) direction is allowed, whereas, for a field polarized along the (1, 0, 0) direction, it is forbidden.In this subsection, we compare the microscopic optical response of GaAs to the excitation by driving fields polarized along the (1, 1, 1) direction and along the (1, 0, 0) direction.
In the presented calculations, we did not take spin-orbit coupling (SOC) into account.Using a smaller Monkhorst-Pack grid (4 Â 4 Â 4), we studied the role of SOC for the calculation of laserdressed GaAs, since SOC plays a role for the electronic structure of GaAs. 66,67In the presence of SOC, bands cannot be assumed to be doubly occupied and we took 120 bands into account.We do observe that spin-up and spin-down components of optically induced charge distributions behave differently, but this effect is small.GaAs is invariant under time-reversal symmetry.Therefore, spin-averaged densities have either only nonzero real or only nonzero imaginary part (see Sec. II C).Comparing spin-averaged charge distributions obtained with and without SOC for driving fields polarized along the (1, 0, 0) direction, we saw that SOC does not play an essential for spinaveraged charge distributions.
The first and second columns of Fig. 2 show, respectively, the first-and second-order oscillations of the electronic state of GaAs driven by an optical field polarized along the (1, 1, 1) direction.The first-order oscillations of laser-driven GaAs comprise the oscillations of the electron current density as Àj 1 ðrÞ cosðxtÞ and the oscillations of the electron density as À. 1 sinðxtÞ.
Figure 2(a) shows the first-order electron current density at xt ¼ 0. The vector field Àj 1 ðrÞ clearly points along the driving-field polarization direction in agreement with the selection rule that the macroscopic first-order polarization of GaAs is aligned with the electric field. 55 The second column of Fig. 2 shows the second-order oscillations of the electronic state that comprise the oscillations of the electron density as . 2 cosð2xtÞ and of the electron current density as Àj 2 sinð2xtÞ.According to the second-order susceptibility tensor of GaAs, 55 the second-order macroscopic polarization driven by an electric field polarized along the (1, 1, 1) direction is also aligned along (1, 1, 1). Figure 2(e) shows the second-order electron density at xt ¼ 0. It also displays a threefold rotational symmetry with respect to the driving-field polarization direction (1, 1, 1).The positive charge alternates with the negative charge along the (1, 1, 1) direction in agreement with the macroscopic polarization aligned along (1, 1, 1).
The magnitude of the second-order electron current density reaches the maximum at xt ¼ p=4 and is shown in Fig. 2(f).It has a very complex structure that is difficult to characterize.We calculate the volume integral of Àj 2 and find that it indeed points in the (1, 1, 1) direction in agreement with the selection rule for the second-order macroscopic polarization.
The first and second columns of Fig. 3 show, respectively, the first-and second-order oscillations of the electronic state of GaAs driven by a field polarized along the (1, 0, 0) direction.Figure 3(a) shows the first-order electron current density at xt ¼ 0, when its magnitude is at a maximum.The vector field clearly points along the driving-field polarization direction (1, 0, 0).This is in agreement with the alignment of the first-order macroscopic polarization of GaAs with the electric field. 55gure 3(b) shows the first-order electron density at xt ¼ p=2.It has twofold rotational symmetry with respect to the direction of the driving-field polarization (1, 0, 0).It is not obvious how the charge distribution in Fig. 3(b) results in the first-order macroscopic polarization along (1, 0, 0), but one may notice that negative charges alter with positive charges in the x direction, when looking at charges around the bottom Ga atoms.
According to the second-order susceptibility tensor of GaAs, 55 its second-order macroscopic polarization is zero for a driving field polarized along the x direction.We obtain that the second-order microscopic optical response is, in fact, nonzero.Figure 3(e) shows the second-order electron density at xt ¼ 0. It also has twofold rotational symmetry about the x axis as . 1 ðrÞ in Fig. 3(b).Figure 3(f) shows the second-order electron current density at xt ¼ p=4.Despite its complex structure, we find in our calculations that its volume integral is indeed zero in agreement with the zero second-order macroscopic polarization.The magnitudes of the second-order electron current density in Figs.3(f) and 3(h) are similar to the magnitudes of the second-order electron current density induced by the field polarized along the (1, 1, 1) direction in Figs.2(f) and 2(h).However, in the first case, the volume integral over the second-order electron current density is negligible and, in the second case, it is considerable.This demonstrates that microscopic optical response resulting in vanishing macroscopic optical response and microscopic optical response resulting in considerable macroscopic optical response can have comparable magnitudes.

IV. CONCLUSIONS
Macroscopic polarization and a HHG spectrum of a crystal are determined by the dipole moment of the electron density.On the atomic scale, optically induced charge distributions go far beyond the concept of a dipole and have a complex spatial structure.This structure has several interesting properties determined by the symmetry of the crystal.Time-reversal symmetry determines the phase of lth-order oscillations of the optically induced charge distribution and the electron current density.We found that components of the electron density oscillate either in phase with the electric field or in phase with the vector potential depending on the parity of the oscillation order.Evenorder charge distributions evolve as harmonics in phase with the vector potential of the optical field, whereas odd-order charge distributions evolve as harmonics in phase with the electric field.For the evenand odd-order amplitudes of the electron current density, this is opposite.Spatial inversion symmetry of the crystal leads to the spatial inversion symmetry of lth-order optically induced charge distributions and electron-current-density amplitudes.Thereby, even-order electron-current-density amplitudes are antisymmetric with respect to the transformation r !Àr, and odd-order electron-current-density amplitudes are symmetric.Even-order optically induced charge distributions are symmetric, and odd-order distributions are antisymmetric.As a result, odd-order optically induced charge distributions are aligned in such a way that macroscopic polarization is induced, and even-order distributions lead to a vanishing macroscopic polarization and, thus, vanishing even-order harmonics.Thus, even when macroscopic optical response is forbidden, charges still rearrange within the unit cell of the crystal.Our analysis within the Floquet-Bloch formalism and calculations in combination with the DFT provide the microscopic insight into optically induced charge distributions due to strong-field optical excitation and provide a deeper understanding of relevant properties of strong-field response such as HHG in crystals.

FIG. 1 .
FIG. 1.The lth-order microscopic optical response of a MgO crystal at different phases of the driving electromagnetic field polarized along the z direction, where l ¼ 1, 2, 3, and 4. A cut of a unit cell centered around the Mg atom is shown.The first column shows the oscillations of the electron density and the electron current density with frequency x at phases xt equal to (a) 0; (b) p=2; (c) p; and (d) 3p=2.Second column corresponds to 2x oscillations at phases xt equal to (e) 0; (f) p=4; (g) p=2; and (h) 3p=4.Third column shows 3x oscillations at phases xt equal to (i) 0; (j) p=6 (k) p=3; and (l) 3p=6.Fourth column shows 4x oscillations with at phases xt equal to (m) 0; (n) p=8; (o) p=4; and (p) 3p=8.The yellow and blue colors represent negative and positive charges, respectively.
Figure 2(b) shows the first-order electron density at xt ¼ 0. It has a threefold rotational symmetry with respect to the driving-field polarization direction (1, 1, 1).The positive charge alternates with the negative charge along the (1, 1, 1) direction.

FIG. 2 .
FIG. 2. The first-and second-order microscopic optical response of a GaAs crystal at different phases of the driving electromagnetic field polarized along the (1, 1, 1) direction.A cut of a unit cell of GaAs centered at the As atom is shown.The first column shows the oscillations of the electron density and the electron current density with the frequency x at phases xt equal to (a) 0; (b) p=2; (c) p; and (d) 3p=2.Second column corresponds to 2x-oscillations at phases xt equal to (e) 0; (f) p=4; (g) p=2; and (h) 3p=4.The yellow and blue colors represent negative and positive charges, respectively.